#ifndef TGraphAnalysis_h
#define TGraphAnalysis_h

#include "TGraph.h"
#include "TMatrixD.h"
#include "RooChi2Var.h
#include "RooMinuit.h
#include "TString.h"
#include "TLine.h"
using namespace RooFit ;


///usage:
///usage: read TGraph *gr to find 
///usage: 	* minimal element
///usage:    	* elements el1,el2 for which Delta1=gr[el1] - gr_minimal=1; Delta2=gr[el2] - gr_minimal=1;
///usage:
///usage:	el1 and el2 might be used to estimate errors (uncorrelated)
///usage:

#include <math.h>
#define PI 3.14159265


void TGraphAnalysis(TGraph* gr,Double_t &el1,Double_t &el2,Double_t &yel1, Double_t &yel2,Double_t & xmin,Double_t & ymin ) ///usage
{
el1=-1e10;
el2=-1e10;
yel1=-1e10;
yel2=-1e10;
xmin=-1e10;
ymin=-1e10;

if (!gr) return;





Double_t _x ;
Double_t _y ;


Double_t _xmin;
Double_t _ymin=1e30;
Int_t _nLess1=0;
Double_t _aLess1[500];
Double_t _bLess1[500];

for (int i=0;i<gr->GetN();i++)
{
gr->GetPoint(i,_x,_y);
if (_y<_ymin) { _ymin=_y;_xmin=_x ;}
}


cout<<"Min 2 "<<_ymin<<"\n";

for (int i=0;i<gr->GetN();i++)
{
gr->GetPoint(i,_x,_y);

//cout<<"_y = "<<_y<<"\n";


if ( (_y-_ymin) <= 1e0 ) {
_aLess1[_nLess1]=_x;
_bLess1[_nLess1]=_y;
_nLess1++;

}

}

cout<<"_nLess1 = "<<_nLess1<<"\n";
cout<<"_xmin = "<<_xmin<<"\n";
cout<<"_ymin = "<<_ymin<<"\n";
cout<<"_aLess1[0] = "<<_aLess1[0]<<"\n";
cout<<"_aLess1[_nLess1-1] = "<<_aLess1[_nLess1-1]<<"\n";

xmin=_xmin;
ymin=_ymin;

if (_aLess1[0]<=_xmin) {el1=_aLess1[0]; yel1=_bLess1[0];}
if (_aLess1[_nLess1-1] >=_xmin) {el2=_aLess1[_nLess1-1];  yel2=_bLess1[_nLess1-1];}



return;

}

void TGraphGetParabola(TGraph* gr, TF1 * & fun) ///usage: 
{
fun=0x0;


if (!gr) return;

Double_t _x ;
Double_t _y ;


Double_t _xmin;
Double_t _ymin=1e30;
Int_t _nLess1=0;
Double_t _aLess1[500];
Double_t _bLess1[500];

for (int i=0;i<gr->GetN();i++)
{
gr->GetPoint(i,_x,_y);
if (_y<_ymin) { _ymin=_y;_xmin=_x ;}
}


for (int i=0;i<gr->GetN();i++)
{
gr->GetPoint(i,_x,_y);

//cout<<"_y = "<<_y<<"\n";


if ( (_y-_ymin) <= 1e2 && _nLess1<499 ) {
_aLess1[_nLess1]=_x;
_bLess1[_nLess1]=_y;
_nLess1++;

}


}

Int_t _tmp[2];
Int_t nd=0;
for (;;)
{
Int_t nd1 =  (Int_t) gRandom->Uniform(0, _nLess1);
if (_aLess1[nd1] != _xmin) { _tmp[nd]=nd1; nd++;}
if (nd==2) break;
}

Double_t el1,el2,el3;
Double_t yel1,yel2,yel3;

el3=_xmin;
yel3=_ymin;

el1=_aLess1[_tmp[0]];
yel1=_bLess1[_tmp[0]];
el2=_aLess1[_tmp[1]];
yel2=_bLess1[_tmp[1]];





///Solve for parabola equaton

TMatrixD m2(3,3);
m2(0,0)=el1*el1;
m2(0,1)=el1;
m2(0,2)=1;

m2(1,0)=el2*el2;
m2(1,1)=el2;
m2(1,2)=1;

m2(2,0)=el3*el3;
m2(2,1)=el3;
m2(2,2)=1;


TVectorD r2(3);

r2(0)=yel1;
r2(1)=yel2;
r2(2)=yel3;

TDecompLU lu2(m2);
lu2.Solve(r2);

Double_t _x1;
Double_t _x2;

gr->GetPoint(0,_x1,_y);
gr->GetPoint(gr->GetN()-1,_x2,_y);



fun = new TF1 ("fun","[0]*x*x+[1]*x+[2]",_x1,_x2);
fun->SetParameter(0,r2(0));
fun->SetParameter(1,r2(1));
fun->SetParameter(2,r2(2));






return;
}


void TEllipseAnalysis(TGraph* gr, Double_t & a, Double_t &b,Double_t xmin,Double_t ymin)
{

/// a-- major axis, b -- minor axis

a=-1e10;
b=-1e10;

if (!gr) return;

Double_t _x ;
Double_t _y ;

///Direct calculation of ellipse's major axis and angle of the rotation 
/// it's valid if we have complete shape
Int_t pos_lft;
Int_t pos_rght;

Double_t _xlft=1e30;
Double_t _ylft=1e30;
Double_t _xrght=-1e30;
Double_t _yrght=-1e30;





for (int i=0;i<gr->GetN();i++)
{
gr->GetPoint(i,_x,_y);
if (_x<_xlft) {pos_lft=i; _ylft=_y;_xlft=_x ;}
if (_x>_xrght) {pos_rght=i; _yrght=_y;_xrght=_x ;}
}

cout<<"_xlft = "<<_xlft<<"\n";
cout<<"_ylft = "<<_ylft<<"\n";
cout<<"_xrght = "<<_xrght<<"\n";
cout<<"_yrght = "<<_yrght<<"\n";

//a=sqrt(pow(_yrght-_ylft,2)+pow(_xrght-_xlft,2));
a=sqrt((_yrght-_ylft)*(_yrght-_ylft)+(_xrght-_xlft)*(_xrght-_xlft));

cout<<"a = "<<a<<"\n";


/*
///Solve for linear equation
/// to find angle of rotation

TMatrixD m2(2,2);

m2(0,0)=_xlft;
m2(0,1)=1;

m2(1,0)=_xrght;
m2(1,1)=1;


TVectorD r2(2);

r2(0)=_ylft;
r2(1)=_yrght;

TDecompLU lu2(m2);
lu2.Solve(r2);


cout<<"r2(0) = "<<r2(0)<<"\n";
cout<<"r2(1) = "<<r2(1)<<"\n";



cout<<"angle = "<<atan(r2(0))*180/PI<<"\n";
cout<<"angle2 = "<<atan2(_yrght-_ylft,_xrght-_xlft)*180/PI<<"\n";

/// The method to find major and minor axis, angle and sigmaX and sigmaY of correlated errors

///Solve for ellipse 
///A*x^2 + B*y^2 +C*x*y + D*x + F*y +G=0; 
/// 6 points!

/// actually it calculates A'*x^2 + B'*y^2 +C'*x*y + D'*x + F'*y =1; A'=A/(-G) etc G is unsovled. It's calculated from knowing the center of ellipes
/// Then A' and B' give a' and b', a=a'/sqrt(G), b=b'/sqrt(G)

/// look at the presentations to understand formulae

TMatrixD m3(5,5);
TVectorD r3(5);

//for (int i=0;i<6;i++)
//r3(i)=0;


for (int i=0;i<5;i++)
{
gr->GetPoint(i,_x,_y);
m3(i,0)=_x*_x;
m3(i,1)=_y*_y;
m3(i,2)=_x*_y;
m3(i,3)=_x;
//m3(i,4)=1;
//r3(i)=_y;

m3(i,4)=_y;
r3(i)=1;
}




//TDecompLU lu3(m3);
TDecompSVD lu3(m3);
lu3.Solve(r3);

///Check consistency on anlge:
Double_t A=r3(0);
Double_t B=r3(1);
Double_t C=r3(2);

Double_t D=0;
if ( A!=B ) D=C/(A-B);

//cout<<"D="<<D<<"\n";

Double_t angle3=asin(D/sqrt(1+D*D))/2;

cout<<"angle3 = "<<asin(D/sqrt(1+D*D))*180/PI/2<<"\n";

Double_t sin2=0;
sin2=D/sqrt(1+D*D);

//cout<<"sin2 = "<<sin2<<"\n";

Double_t aInverse2=0;
Double_t bInverse2=0;

///Dependence on scale factor F!!! (or G )
if (sin2!=0) aInverse2=fabs(C/sin2+A+B)/2.;
if (sin2!=0) bInverse2=fabs(-C/sin2+A+B)/2.;

cout<<"aInverse2 = "<<aInverse2<<"\n";
cout<<"bInverse2 = "<<bInverse2<<"\n";


if (aInverse2!=0) cout<<"a' = "<<sqrt(1e0/aInverse2)<<"\n";
if (bInverse2!=0) cout<<"b' = "<<sqrt(1e0/bInverse2)<<"\n";



cout<<"G' = "<<A*xmin*xmin+B*ymin*ymin+C*xmin*ymin-1<<"\n";
cout<<"G = "<<-1e0/(A*xmin*xmin+B*ymin*ymin+C*xmin*ymin+1)<<"\n";

Double_t G=-1e0/(A*xmin*xmin+B*ymin*ymin+C*xmin*ymin+1);


if (aInverse2!=0 && G!=0) cout<<"a = "<<sqrt(1e0/aInverse2/fabs(G))<<"\n";
if (bInverse2!=0 && G!=0) cout<<"b = "<<sqrt(1e0/bInverse2/fabs(G))<<"\n";




//a=sqrt(1e0/aInverse2);
//a=sqrt(1e0/bInverse2);
a=sqrt(1e0/aInverse2/fabs(G));
b=sqrt(1e0/bInverse2/fabs(G));

Double_t SigmaX;
Double_t SigmaY;






SigmaX=sqrt(sin(angle3)*sin(angle3)*b*b+cos(angle3)*cos(angle3)*a*a);
SigmaY=sqrt(sin(angle3)*sin(angle3)*a*a+cos(angle3)*cos(angle3)*b*b);

 cout<<"SigmaX = "<<SigmaX<<"\n";
 cout<<"SigmaY = "<<SigmaY<<"\n";

*/

///The direct scanning of the distance from center to each point. Then major/minor and angle can be approximated

Double_t _xa;
Double_t _ya;
Double_t _xb;
Double_t _yb;
Double_t _adist=-1e30;
Double_t _bdist=1e30;



for (int i=0;i<gr->GetN();i++)
{
gr->GetPoint(i,_x,_y);

///here some cleaning of result
if (_x<1 || _y<1) continue;

Double_t _distance =sqrt((xmin-_x)*(xmin-_x)+(ymin-_y)*(ymin-_y));
if (_adist<_distance) { _adist=_distance;_xa=_x, _ya=_y;}
if (_bdist>_distance) { _bdist=_distance;_xb=_x, _yb=_y;}
}

if (xmin<1) _adist=_xa;
if (ymin<1) _bdist=_yb;


cout<<"tan = "<<(_ya-ymin)/(_xa-xmin)<<"\n";

Double_t angle4 = atan2(_ya-ymin,_xa-xmin);

cout<<"angle4 = "<<atan2(_ya-ymin,_xa-xmin)*180/PI<<"\n";

cout<<"a = "<<_adist<<"\n";
cout<<"b = "<<_bdist<<"\n";

Double_t SigmaX_2=sqrt(sin(angle4)*sin(angle4)*_bdist*_bdist+cos(angle4)*cos(angle4)*_adist*_adist);
Double_t SigmaY_2=sqrt(sin(angle4)*sin(angle4)*_adist*_adist+cos(angle4)*cos(angle4)*_bdist*_bdist);

 cout<<"SigmaX = "<<SigmaX_2<<"\n";
 cout<<"SigmaY = "<<SigmaY_2<<"\n";






return; 
}

TGraph * TMinuitAnalysis(TMinuit* ptMinuit,Int_t paramN,Double_t elfrom,Double_t elto,Double_t min,Double_t max,Double_t &el1,Double_t &el2,Double_t &yel1, Double_t &yel2,Double_t & xmin,Double_t & ymin ) ///usage
{

const Int_t _nStep=300;
el1=-1e10;
el2=-1e10;
yel1=-1e10;
yel2=-1e10;
xmin=-1e10;
ymin=-1e10;

Double_t _Chi2Vals[_nStep];
Double_t _ParVals[_nStep];


Double_t arglist[10];
Int_t ierflg = 0;


if (!ptMinuit) return;

///Supprese printings
ptMinuit->mnexcm("SET NOWarnings", arglist ,0,ierflg);
ptMinuit->SetPrintLevel(-1);


///Store initial values
Int_t npars = ptMinuit->GetNumPars();

Double_t _vals0[100];
Double_t _errs[100];
Double_t _vals1[100];

for (int i=0;i<npars;i++)
	 ptMinuit->GetParameter(i, _vals0[i],_errs[i]);

Double_t _fval;
ptMinuit->Eval(npars, 0, _fval, _vals0,ierflg);
cout<<"Current Chi2 is "<<_fval<<"\n";

if (elfrom<min) elfrom=min;
if (elto>max) elto=max;


Double_t _step=(elto - elfrom)/_nStep;

arglist[0] =  paramN+1;
ptMinuit->mnexcm("FIX", arglist ,1,ierflg);

Int_t index=0;

///performe stepping 
for (double x=elfrom; x<elto+_step; x+=_step)
{

 arglist[0] =  paramN+1;
 arglist[1] =  x;
 ptMinuit->mnexcm("SET PAR", arglist ,2,ierflg);

  arglist[0] = 3000;
  arglist[1] = 1; ///tolerance
  ptMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);

for (int i=0;i<npars;i++)
         ptMinuit->GetParameter(i, _vals1[i],_errs[i]);
  ptMinuit->Eval(npars, 0, _fval, _vals1,ierflg);

//cout<<"Current _fval = "<<_fval<<"\n";

if (index<_nStep) {
 _Chi2Vals[index]= _fval ;
 _ParVals[index]=x;

}

index++;

}


gr = new TGraph(_nStep,&_ParVals[0],&_Chi2Vals[0]);
TGraphAnalysis( gr,el1,el2,yel1,yel2, xmin, ymin );


///Restore all initial vals and clean memory
arglist[0] =  paramN+1;
ptMinuit->mnexcm("RELease", arglist ,1,ierflg);

for (int i=0;i<npars;i++)
{
 arglist[0] =  i+1;
 arglist[1] =  _vals0[i];
 ptMinuit->mnexcm("SET PAR", arglist ,2,ierflg);

}

 /// Restore HESSE
  arglist[0] = 3000;
  arglist[1] = 1; ///tolerance
  ptMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);

//delete gr;


//cout<<"gr = "<<gr<<endl;


return gr;
}

TGraph * TMinuitCreateContour(TMinuit* ptMinuit,Int_t paramN1,Int_t paramN2,TList *arrayOfNames,Double_t scanRange,Double_t sigmas,Double_t min,Double_t max ) ///usage:
{

const Int_t nRun=10000;

Double_t arglist[10];
Int_t ierflg = 0;


if (!(ptMinuit && arrayOfNames)) return;

TString name1 = arrayOfNames->At(paramN1+1)->GetName();
TString name2 = arrayOfNames->At(paramN2+1)->GetName();

///Suppress printings

ptMinuit->mnexcm("SET NOWarnings", arglist ,0,ierflg);
ptMinuit->SetPrintLevel(-1);


///Store initial values
Int_t npars = ptMinuit->GetNumPars();

Double_t _vals0[40];
Double_t _errs0[40];
Double_t _errs1[40];
Double_t _vals1[40];

for (int i=0;i<npars;i++)
	 ptMinuit->GetParameter(i, _vals0[i],_errs0[i]);


for (int ii=0;ii<npars;ii++){

        _vals1[ii]=_vals0[ii];
        _errs1[ii]=_errs0[ii];

}


Double_t _fvalMin;
Double_t _fval;
ptMinuit->Eval(npars, 0, _fvalMin, _vals0,ierflg);

cout<<"Current Min Chi2 is "<<_fvalMin<<"\n";


const Int_t nBin2D=30;
Double_t minX=0;
Double_t maxX=0;
Double_t minY=0;
Double_t maxY=0;


minX=_vals1[paramN1]-scanRange*_errs1[paramN1];
maxX=_vals1[paramN1]+scanRange*_errs1[paramN1];


if (minX<min) minX=min;
if (maxX>max) maxX=max;

minY=_vals1[paramN2]-scanRange*_errs1[paramN2];
maxY=_vals1[paramN2]+scanRange*_errs1[paramN2];


if (minY<min) minY=min;
if (maxY>max) maxY=max;

TH2D * Hist2D = new TH2D(Form("H2CONT_%s_%s",name1.Data(),name2.Data()),Form("H2CONT_%s_%s",name1.Data(),name2.Data()),nBin2D,minX,maxX,nBin2D,minY,maxY);
Hist2D->SetDirectory(gROOT);

arglist[0] =  paramN1+1;
ptMinuit->mnexcm("FIX", arglist ,1,ierflg);
arglist[0] =  paramN2+1;
ptMinuit->mnexcm("FIX", arglist ,1,ierflg);

for (int k=0;k<nRun;k++)
{

Double_t _X = gRandom->Uniform(minX,maxX);
Double_t _Y = gRandom->Uniform(minY,maxY);


 arglist[0] =  paramN1+1;
 arglist[1] =  _X;
 ptMinuit->mnexcm("SET PAR", arglist ,2,ierflg);
 arglist[0] =  paramN2+1;
 arglist[1] =  _Y;
 ptMinuit->mnexcm("SET PAR", arglist ,2,ierflg);


 arglist[0] = 500;
  arglist[1] = 1; ///tolerance
  ptMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);

for (int kk=0;kk<npars;kk++)
        ptMinuit->GetParameter(kk, _vals0[kk],_errs0[kk]);

ptMinuit->Eval(npars, 0, _fval, _vals0,ierflg);
for (int iii=0;iii<npars;iii++)
ptMinuit->DefineParameter(iii,Form("yield_%s",arrayOfNames->At(iii+1)->GetName()),_vals1[iii],_errs1[iii],min,max);

if (_fval<=_fvalMin+sigmas*sigmas) {
 Hist2D->Fill(_X,_Y);
}


} /// for nRun


///Restore all initial vals
arglist[0] =   paramN1+1;
ptMinuit->mnexcm("RELease", arglist ,1,ierflg);
arglist[0] =   paramN2+1;
ptMinuit->mnexcm("RELease", arglist ,1,ierflg);


Hist2D->Draw("SCAT");

///Calculate TGraph for analysis and plottings

Double_t DownLimitX[nBin2D];
Double_t DownLimitY[nBin2D];
Double_t UpLimitX[nBin2D];
Double_t UpLimitY[nBin2D];
Int_t numPointsDown=0;
Int_t numPointsUp=0;
Int_t row=0;

for (int k1=1;k1<=nBin2D;k1++)
{
Int_t numInternalPoints=0;
for (int k2=1;k2<=nBin2D;k2++)
{

Double_t _cont = Hist2D->GetBinContent(k1,k2);
Double_t xcnt;
Double_t ycnt;

if (_cont!=0) {
numInternalPoints++;
TAxis *xaxis = Hist2D->GetXaxis();
TAxis *yaxis = Hist2D->GetYaxis();
xcnt = xaxis->GetBinCenter(k1);
ycnt = yaxis->GetBinCenter(k2);

if (numInternalPoints==1)
{
        DownLimitX[numPointsDown]=xcnt;
        DownLimitY[numPointsDown]=ycnt;
        numPointsDown++;
}

if (numInternalPoints>1)
{
        UpLimitX[numPointsUp]=xcnt;
        UpLimitY[numPointsUp]=ycnt;

}

if (k2==nBin2D && numInternalPoints>1) numPointsUp++;

} /// for k2
} /// for k1


TGraph *gr3 =0x0;
if (numPointsDown+numPointsUp>0) {
gr3= new TGraph(numPointsDown+numPointsUp);
for (int k1=0; k1<numPointsDown;k1++)
gr3->SetPoint(k1,DownLimitX[k1],DownLimitY[k1]);
for (int k2=0;k2<numPointsUp;k2++)
gr3->SetPoint(numPointsDown+k2,UpLimitX[k2],UpLimitY[k2]);

gr3->SetName(Form("gr3_CONT_%s_%s",name1.Data(),name2.Data()));
gr3->SetTitle(Form("CONT_%s_%s",name1.Data(),name2.Data()));
}

///Allow printings
ptMinuit->mnexcm("SET WARnings", arglist ,0,ierflg);
ptMinuit->SetPrintLevel();


for (int ii=0;ii<npars;ii++)
{
ptMinuit->DefineParameter(ii,Form("yield_%s",arrayOfNames->At(ii+1)->GetName()),_vals1[ii],_errs1[ii],min,max);

}



 /// Restore HESSE
  arglist[0] = 3000;
  arglist[1] = 1; ///tolerance
  ptMinuit->mnexcm("MINOS", arglist ,2,ierflg);
}

return gr3;
}



///usage: Plotting functions

void TGraphLinePlot(TGraph * gr,Color_t col=kBlue,Double_t el1,Double_t el2,Double_t xmin,TAxis * yAx=0 ,TString options="",Double_t offset=0.0) ///usage: plot lines with parabola 
{
	if (!gr) return;

if (options.Length()>0) gr->Draw(options);
else gr->Draw("a*");


double _m1; 
double _m2;

if ( !(yAx)) {
_m1= gr->GetYaxis()->GetXmin();
_m2 = gr->GetYaxis()->GetXmax();
}
else
{
_m1= yAx->GetXmin();
_m2= yAx->GetXmax();

}

TLine *l1 = new TLine(el1,_m1,el1,_m2);
TLine *l2 = new TLine(el2,_m1,el2,_m2);
TLine *l3 = new TLine(xmin,_m1,xmin,_m2);

l1->SetLineColor(col);
l2->SetLineColor(col);
l3->SetLineColor(kGreen);

l1->Draw("same");
l2->Draw("same");
l3->Draw("same");

/*
l1->SetBit(TLine::kLineNDC);
Double_t x1=l1->GetX1();

cout<<"X1 "<<l1->GetX1();
cout<<"X2 "<<l1->GetX2();
cout<<"Y1 "<<l1->GetY1();
cout<<"Y2 "<<l1->GetY2();
*/


 TText *label1 = new TText();
 label1 -> SetTextFont(1);
 label1 -> SetTextColor(col);
 label1 -> SetTextSize(0.02);
 label1 -> SetTextAlign(22);
 label1 -> SetTextAngle(0);
if (el1==0) el1=0.1;
if (offset==0.0)
 label1 -> DrawText(el1,-20*_m1 ,Form(" %.1f",el1));
else
 label1 -> DrawText(el1,_m1*(1.+offset) ,Form(" %.1f",el1));

 TText *label2 = new TText();
 label2 -> SetTextFont(1);
 label2 -> SetTextColor(col);
 label2 -> SetTextSize(0.02);
 label2 -> SetTextAlign(22);
 label2 -> SetTextAngle(0);
if (el2==0) el2=0.1;

if (offset==0.0)
 label2 -> DrawText(el2,-30*_m1 ,Form(" %.1f",el2));
else
label2 -> DrawText(el2,_m1*(1.+offset*2) ,Form(" %.1f",el2));


 TText *label3 = new TText();
 label3 -> SetTextFont(1);
 label3 -> SetTextColor(kGreen);
 label3 -> SetTextSize(0.02);
 label3 -> SetTextAlign(22);
 label3 -> SetTextAngle(0);
if (xmin==0) xmin=0.1;

if (offset==0.0)
 label3 -> DrawText(xmin,-25*_m1 ,Form(" %.1f",xmin));
else
label3 -> DrawText(xmin,_m1*(1.+offset*3) ,Form(" %.1f",xmin));



}

void TGraphEllipsePlot(TGraph * gr,Color_t col1=kBlue,Color_t col2=kRed,Double_t xmin,Double_t ymin,TString options="") ///usage: plot lines with parabola 
{
	if (!gr) return;

if (options.Length()>0) gr->Draw(options);
else gr->Draw("a*");

///The direct scanning of the distance from center to each point. Then major/minor and angle can be approximated

Double_t _xa,_x;
Double_t _ya,_y;
Double_t _xb;
Double_t _yb;
Double_t _adist=-1e30;
Double_t _bdist=1e30;

Double_t a,b;

TEllipseAnalysis(gr,a,b,xmin,ymin);




for (int i=0;i<gr->GetN();i++)
{
gr->GetPoint(i,_x,_y);

///here some cleaning of result
if (_x<1 || _y<1) continue;

Double_t _distance =sqrt((xmin-_x)*(xmin-_x)+(ymin-_y)*(ymin-_y));
if (_adist<_distance) { _adist=_distance;_xa=_x, _ya=_y;}
if (_bdist>_distance) { _bdist=_distance;_xb=_x, _yb=_y;

//cout<<"\n"<<"_bdist = "<<_bdist<<"\n";
//cout<<"\n"<<"_xb = "<<_xb<<"\n";
//cout<<"\n"<<"_yb = "<<_yb<<"\n";
}
}

if (xmin<1) _adist=_xa;
if (ymin<1) _bdist=_yb;


TLine *l1 = new TLine(_xa,_ya,xmin,ymin);
TLine *l2 = new TLine(_xb,_yb,xmin,ymin);

l1->SetLineColor(col1);
l2->SetLineColor(col2);

l1->Draw("same");
l2->Draw("same");
TMarker *m = new TMarker(xmin, ymin, 29);
m->Draw("same");


 TText *label1 = new TText();
 label1 -> SetTextFont(1);
 label1 -> SetTextColor(col1);
 label1 -> SetTextSize(0.02);
 label1 -> SetTextAlign(22);
 label1 -> SetTextAngle(0);
//if (_xa==0) _xa=0.1;
 label1 -> DrawText(_xa,_ya ,Form(" %.1f",_adist));

 TText *label2 = new TText();
 label2 -> SetTextFont(1);
 label2 -> SetTextColor(col2);
 label2 -> SetTextSize(0.02);
 label2 -> SetTextAlign(22);
 label2 -> SetTextAngle(0);
//if (_xb==0) _xb=0.1;

 label2 -> DrawText(_xb,_yb ,Form(" %.1f",_bdist));



 TText *label3 = new TText();
 label3 -> SetTextFont(1);
 label3 -> SetTextColor(kGreen);
 label3 -> SetTextSize(0.02);
 label3 -> SetTextAlign(22);
 label3 -> SetTextAngle(0);
//if (xmin==0) xmin=0.1;
//if (ymin==0) ymin=0.1;
 label3 -> DrawText(xmin,ymin ,Form(" %.1f",xmin));


return;
}




///usage: RooFit stuff 
TGraph * RooFitChi2Analysis(RooChi2Var* chi2, RooRealVar* varScan,Double_t sigmas,Double_t min,Double_t max,Double_t &el1,Double_t &el2,Double_t &yel1, Double_t &yel2,Double_t & xmin,Double_t & ymin ) ///usage
{


TGraph * gr8=0x0;

if ( !(chi2 && varScan) ) return gr8;

Double_t _val =varScan->getVal();
Double_t _err =varScan->getError();


const Int_t _nStep=300;
Double_t _Chi2Vals[_nStep];
Double_t _ParVals[_nStep];
Int_t index=0;

Double_t elfrom= varScan->getVal() - sigmas *varScan->getError();
Double_t elto= varScan->getVal() + sigmas *varScan->getError();


if (elfrom<min) elfrom=min;
if (elto>max) elto=max;

Double_t _step=(elto - elfrom)/_nStep;

///performe stepping
for (double x=elfrom; x<elto+_step; x+=_step)
{
varScan-> setVal( x ) ;

if (index<_nStep) {
 _Chi2Vals[index]= chi2->getVal() ;
 _ParVals[index]=x;
}

index++;


}

TGraph * gr8 = new TGraph(_nStep,&_ParVals[0],&_Chi2Vals[0]);

TGraphAnalysis( gr8,el1,el2,yel1,yel2, xmin, ymin );

varScan-> setVal( _val) ;
varScan-> setError( _err) ;


return gr8;
}

///usage:  lstVars contains a data yield as the first element
TGraph * RooFitMinuitChi2Analysis(RooMinuit *minuit,RooChi2Var* chi2, RooRealVar* varScan,TList * lstVars,Double_t sigmas,Double_t min,Double_t max,Double_t &el1,Double_t &el2,Double_t &yel1, Double_t &yel2,Double_t & xmin,Double_t & ymin ) ///usage
{


TGraph * gr8=0x0;

if ( !(minuit && chi2 && varScan && lstVars) ) return gr8;

//Store all vals

Double_t _vals[40];
Double_t _errs[40];

Int_t nPars =  lstVars->GetEntries();


for (int i=1;i<nPars;i++)
{
	_vals[i-1]=((RooRealVar*)arrayOfYields->At(i))->getVal();
	_errs[i-1]=((RooRealVar*)arrayOfYields->At(i))->getError();

}
Int_t printLevel = minuit->setPrintLevel(-1);
Int_t warnLevel = minuit->setWarnLevel(-1);

///Double_t _val =varScan->getVal();


const Int_t _nStep=300;
Double_t _Chi2Vals[_nStep];
Double_t _ParVals[_nStep];
Int_t index=0;

Double_t elfrom= varScan->getVal() - sigmas *varScan->getError();
Double_t elto= varScan->getVal() + sigmas *varScan->getError();


if (elfrom<min) elfrom=min;
if (elto>max) elto=max;

Double_t _step=(elto - elfrom)/_nStep;

varScan-> setConstant( kTRUE ) ;

///performe stepping
for (double x=elfrom; x<elto+_step; x+=_step)
{

///Restore initial values each step
for (int i=1;i<nPars;i++)
{
	if (arrayOfYields->At(i) == varScan) continue;
        ((RooRealVar*)arrayOfYields->At(i))->setVal(_vals[i-1]);
    	((RooRealVar*)arrayOfYields->At(i))->setError( _errs[i-1]);

}



varScan-> setVal( x ) ;

minuit->migrad() ;
minuit->hesse() ;

if (index<_nStep) {
 _Chi2Vals[index]= chi2->getVal() ;
 _ParVals[index]=x;
}

index++;


}

TGraph * gr8 = new TGraph(_nStep,&_ParVals[0],&_Chi2Vals[0]);

TGraphAnalysis( gr8,el1,el2,yel1,yel2, xmin, ymin );


///Restore initial values each step
for (int i=1;i<nPars;i++)
{
        ((RooRealVar*)arrayOfYields->At(i))->setVal(_vals[i-1]);
        ((RooRealVar*)arrayOfYields->At(i))->setError( _errs[i-1]);

}

varScan-> setConstant( kFALSE ) ;

minuit->setPrintLevel(printLevel);
minuit->setWarnLevel(warnLevel);


return gr8;
}


TGraph * RooFitCreateContour(RooMinuit *minuit,RooChi2Var* chi2, RooRealVar* varScan1,  RooRealVar* varScan2,TList * lstVars,Double_t scanRange,Double_t sigmas,Double_t min,Double_t max ) ///usage	
{
TGraph *gr=0x0;

if (!(minuit && chi2 && varScan1 && varScan2 && lstVars)) return gr;

const Int_t nRun=10000;
const Int_t nBin2D=30;
Double_t minX=0;
Double_t maxX=0;
Double_t minY=0;
Double_t maxY=0;

Double_t _vals[40];
Double_t _errs[40];

Int_t nPars =  lstVars->GetEntries();


for (int i=1;i<nPars;i++)
{
        _vals[i-1]=((RooRealVar*)arrayOfYields->At(i))->getVal();
        _errs[i-1]=((RooRealVar*)arrayOfYields->At(i))->getError();

	

}

Int_t printLevel = minuit->setPrintLevel(-1);
Int_t warnLevel = minuit->setWarnLevel(-1);




Double_t _fvalMin = chi2->getVal();
minX =  varScan1->getVal() - scanRange *varScan1->getError();
maxX= varScan1->getVal() + scanRange *varScan1->getError();



if (minX<min) minX=min;
if (maxX>max) maxX=max;

minY =  varScan2->getVal() - scanRange *varScan2->getError();
maxY= varScan2->getVal() + scanRange *varScan2->getError();



if (minY<min) minY=min;
if (maxY>max) maxY=max;

varScan1-> setConstant( kTRUE ) ;
varScan2-> setConstant( kTRUE ) ;

TH2D * Hist2D = new TH2D(Form("H2CONT_%s_%s",varScan1->GetName(),varScan2->GetName()),Form("H2CONT_%s_%s",varScan1->GetName(),varScan2->GetName()),nBin2D,minX,maxX,nBin2D,minY,maxY);
Hist2D->SetDirectory(gROOT);


for (int k=0;k<nRun;k++)
{

Double_t _X = gRandom->Uniform(minX,maxX);
Double_t _Y = gRandom->Uniform(minY,maxY);

///Restore initial values each step
for (int i=1;i<nPars;i++)
{
        if (arrayOfYields->At(i) == varScan1) continue;
        if (arrayOfYields->At(i) == varScan2) continue;
        ((RooRealVar*)arrayOfYields->At(i))->setVal(_vals[i-1]);
        ((RooRealVar*)arrayOfYields->At(i))->setError( _errs[i-1]);

}



varScan1-> setVal( _X ) ;
varScan2-> setVal( _Y ) ;

minuit->migrad() ;
minuit->hesse() ;

Double_t _fval = chi2->getVal();

if (_fval<=_fvalMin +sigmas*sigmas ) {
 Hist2D->Fill(_X,_Y);
}

} //for nRun

Hist2D->Draw("SCAT");

///Calculate TGraph for analysis and plottings

Double_t DownLimitX[nBin2D];
Double_t DownLimitY[nBin2D];
Double_t UpLimitX[nBin2D];
Double_t UpLimitY[nBin2D];
Int_t numPointsDown=0;
Int_t numPointsUp=0;
Int_t row=0;

for (int k1=1;k1<=nBin2D;k1++)
{
Int_t numInternalPoints=0;
for (int k2=1;k2<=nBin2D;k2++)
{


Double_t _cont = Hist2D->GetBinContent(k1,k2);
Double_t xcnt;
Double_t ycnt;

if (_cont!=0) {
numInternalPoints++;
TAxis *xaxis = Hist2D->GetXaxis();
TAxis *yaxis = Hist2D->GetYaxis();
xcnt = xaxis->GetBinCenter(k1);
ycnt = yaxis->GetBinCenter(k2);

if (numInternalPoints==1)
{
        DownLimitX[numPointsDown]=xcnt;
        DownLimitY[numPointsDown]=ycnt;
        numPointsDown++;
}

if (numInternalPoints>1)
{
        UpLimitX[numPointsUp]=xcnt;
        UpLimitY[numPointsUp]=ycnt;

}

} //if _cont!=0

if (k2==nBin2D && numInternalPoints>1) numPointsUp++;

} /// for k2
} /// for k1

if (numPointsDown+numPointsUp>0) {
gr= new TGraph(numPointsDown+numPointsUp);

for (int k1=0; k1<numPointsDown;k1++)
gr->SetPoint(k1,DownLimitX[k1],DownLimitY[k1]);

for (int k2=0;k2<numPointsUp;k2++)
gr->SetPoint(numPointsDown+k2,UpLimitX[k2],UpLimitY[k2]);

gr->SetName(Form("gr3_CONT_%s_%s",varScan1->GetName(), varScan2->GetName()));
gr->SetTitle(Form("CONT_%s_%s",varScan1->GetName(), varScan2->GetName()));
}

///Restore initial values each step
for (int i=1;i<nPars;i++)
{
        ((RooRealVar*)arrayOfYields->At(i))->setVal(_vals[i-1]);
        ((RooRealVar*)arrayOfYields->At(i))->setError( _errs[i-1]);

}

varScan1-> setConstant( kFALSE ) ;
varScan2-> setConstant( kFALSE ) ;

minuit->setPrintLevel(printLevel);
minuit->setWarnLevel(warnLevel);


return gr;
}

#endif
